#R version 4.4.3 (2025-02-28 ucrt)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 11 x64 

library(dplyr) #1.1.4
library(ggplot2) #3.5.1
library(readr) #2.1.5
library(haven) #2.5.4
library(estimatr) #1.0.6
library(broom.mixed) #0.2.9.6
library(lme4) #1.1-37
library(countrycode) #1.6.1
library(forcats) #1.0.0
library(modelsummary) #2.3.0

`%notin%` <- function(x,y) !(x %in% y) 

setwd("C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/cps/replication")

#load in primary data
gwf <- read_csv("data/gwf_panel.csv") |>  
  select(countryname, year, v2xps_party, v2x_polyarchy, gwf_duration, gwf_regime, gwf_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm) |>  
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown")) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.nan(.) | gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.na(.), 0, .)) |>  
  filter(gwf_prior %in% c("Party", "Military", "Personal") & gwf_regime == "democracy") |>  
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) |>  
  mutate(gwf_prior = factor(gwf_prior, levels = c("Party", "Personal", "Military"))) |>  
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) |> 
  group_by(ccode) |> 
  mutate(ldv = lag(v2xps_party), reactive = case_when(ruling == 1 ~ 0, ruling == 0 ~ 1))

#filter out NA for each model
gwf2 <- gwf |> select(v2xps_party, v2x_polyarchy, lgdp, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()
gwf3 <- gwf |> select(countryname, year, v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, ccode) |> na.omit()
gwf4 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, gwf_pduration, ccode) |> na.omit()
gwf4 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, gwf_pduration, regime_elecs, ccode) |> na.omit()
gwf5 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, gwf_number, fsu, ussr_sat, av_dm, ccode) |> na.omit()
gwf6 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, gwf_number, fsu, ussr_sat, av_dm, ccode) |> na.omit()
gwf7 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()
gwf8 <- gwf |> select(ldv, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()
gwf9 <- gwf |> select(year, ccode, ldv, gwf_prior, pi1, prior_pi, reactive, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()

#run models -- PCSE for models 7, 8, 9 are found in Table2.do
m1 <- lmer(v2xps_party ~ + (1 | ccode), data = gwf, REML=FALSE) #
m2 <- lmer(v2xps_party ~ v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf2, REML=FALSE) #
m3 <- lmer(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + (1 | ccode), data = gwf3, REML=FALSE) # 
m4 <- lmer(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + gwf_pduration + (1 | ccode), data = gwf4, REML=FALSE) #
m5 <- lmer(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + av_dm + (1 | ccode), data = gwf5, REML=FALSE)
m7 <- lmer(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf7, REML=FALSE)
m8 <- lm(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = gwf7) 
sd(residuals(m8))
m9 <- lm(ldv ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = gwf8) 
sd(residuals(m9))
m10 <- lm(ldv ~ gwf_prior + pi1 + prior_pi + reactive + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + regime_elecs + av_dm, data = gwf9) 
sd(residuals(m10))

modellist <- list("ML1" = m1, #just baseline
                  "ML2" = m2, #shows basline with time-variant indicators
                  "ML3" = m3, #just time invariant - prior duration 
                  "ML4" = m4, #add in prior duration, not much change 
                  "ML5" = m5, #full minus prior duration & # of elections, prior_pi drops 
                  "ML7" = m7,
                  "LM" = m8,
                  "LDV" = m9,
                  "LDV" = m10)

modelsummary(modellist, stars = TRUE, gof_omit = "R2|AIC|RM|F|Log.Lik", coef_map = c("pi1" = "Incumbent PI", "prior_pi" = "Prior PI", "asp_dummy" = "ASP",
                                                                                     "reactive" = "Reactive ASP",
                                                                                                              "gwf_priorPersonal" = "Personalist", "gwf_priorMilitary" = "Military",
                                                                                                              "gwf_pduration" = "Prior Duration", "regime_elecs" = "Elections", 
                                                                                                              "SD (Intercept ccode)" = "SD: Country", "SD (Observations)" = "SD: Residuals"), output = "latex")